library(haven)
library(lattice)
library(misty)
library(Matrix)
library(lme4)
library(msm)
library(sandwich)
library(sjPlot)
library(backports)     
library(carData)
library(effects)       
library(ggplot2)       
library(interactions) 
library(lme4)          
library(lmerTest)      
library(psych)         
library(plyr) 
library(ggplot2)
library(ggpubr)
library(ggplot2)

#################################################################################################################
#RMSSD
########################################################################################################################

dataset[dataset == 777] <- NA
dataset[dataset == 999] <- NA

#Base model
dataset_rmssdbase <- dataset[ , c("id", "stressmarker", "resting_rmssd", "reactivity_rmssd", "resting_met", "rmssdrecovery_MET", "rmssd_recovery", "conttime")]
dataset_rmssdbase <- dataset_rmssdbase[complete.cases(dataset_rmssdbase),]

#RMSSD reactivity 
dataset_rmssdbase$rmssd_reactivity <- (dataset_rmssdbase$reactivity_rmssd - dataset_rmssdbase$resting_rmssd)

#RMSSD reactivity dichotom
dataset_rmssdbase$rmssd_reactivity_dichotom <- dataset_rmssdbase$rmssd_reactivity
dataset_rmssdbase$rmssd_reactivity_dichotom [which(dataset_rmssdbase$rmssd_reactivity_dichotom < 0)] <- 0
dataset_rmssdbase$rmssd_reactivity_dichotom [which(dataset_rmssdbase$rmssd_reactivity_dichotom > 0)] <- 1

subsetlowreact <- subset(dataset_rmssdbase, rmssd_reactivity_dichotom==0)
subsethighreact <- subset(dataset_rmssdbase, rmssd_reactivity_dichotom>0)

summary(subsetlowreact$rmssd_reactivity)
summary(subsethighreact$rmssd_reactivity)

#RMSSD recovery 
dataset_rmssdbase$rmssd_recoveryno999 <- dataset_rmssdbase$rmssd_recovery
dataset_rmssdbase$rmssd_recoveryno999 [which(dataset_rmssdbase$rmssd_recoveryno999 > 900)] <- NA
dataset_rmssdbase <- dataset_rmssdbase[complete.cases(dataset_rmssdbase),]

#Log transform RMSSD parameters
#RMSSD resting
dataset_rmssdbase$resting_rmssd_log <- log(dataset_rmssdbase$resting_rmssd)

#Center independent variables
#Resting MET
dataset_rmssdbase$resting_met_cent <- scale(dataset_rmssdbase$resting_met, scale = FALSE)
#Recovery MET
dataset_rmssdbase$rmssdrecovery_MET_cent <- scale(dataset_rmssdbase$rmssdrecovery_MET, scale = FALSE)
#RMSSD Resting
dataset_rmssdbase$resting_rmssd_log_cent <- center(dataset_rmssdbase$resting_rmssd_log, type = c("CWC"), cluster = dataset_rmssdbase$id, value = NULL, as.na = NULL,
                                                   check = TRUE)
#Conttime_hour
dataset_rmssdbase$conttime_cent <- scale(dataset_rmssdbase$conttime, scale = FALSE)
dataset_rmssdbase$conttime_cent <- dataset_rmssdbase$conttime_cent/24

#outliers
#RMSSD resting
dataset_rmssdbase$resting_rmssd_log_cent[which(dataset_rmssdbase$resting_rmssd_log_cent < (mean(dataset_rmssdbase$resting_rmssd_log_cent) - 3*sd(dataset_rmssdbase$resting_rmssd_log_cent)))] <- (mean(dataset_rmssdbase$resting_rmssd_log_cent) - 3*sd(dataset_rmssdbase$resting_rmssd_log_cent))
dataset_rmssdbase$resting_rmssd_log_cent[which(dataset_rmssdbase$resting_rmssd_log_cent > (mean(dataset_rmssdbase$resting_rmssd_log_cent) + 3*sd(dataset_rmssdbase$resting_rmssd_log_cent)))] <- (mean(dataset_rmssdbase$resting_rmssd_log_cent) + 3*sd(dataset_rmssdbase$resting_rmssd_log_cent))
#RMSSD MET
dataset_rmssdbase$resting_met_cent[which(dataset_rmssdbase$resting_met_cent < (mean(dataset_rmssdbase$resting_met_cent) - 3*sd(dataset_rmssdbase$resting_met_cent)))] <- (mean(dataset_rmssdbase$resting_met_cent) - 3*sd(dataset_rmssdbase$resting_met_cent))
dataset_rmssdbase$resting_met_cent[which(dataset_rmssdbase$resting_met_cent > (mean(dataset_rmssdbase$resting_met_cent) + 3*sd(dataset_rmssdbase$resting_met_cent)))] <- (mean(dataset_rmssdbase$resting_met_cent) + 3*sd(dataset_rmssdbase$resting_met_cent))
#RMSSD recovery MET
dataset_rmssdbase$rmssdrecovery_MET_cent[which(dataset_rmssdbase$rmssdrecovery_MET_cent < (mean(dataset_rmssdbase$rmssdrecovery_MET_cent) - 3*sd(dataset_rmssdbase$rmssdrecovery_MET_cent)))] <- (mean(dataset_rmssdbase$rmssdrecovery_MET_cent) - 3*sd(dataset_rmssdbase$rmssdrecovery_MET_cent))
dataset_rmssdbase$rmssdrecovery_MET_cent[which(dataset_rmssdbase$rmssdrecovery_MET_cent > (mean(dataset_rmssdbase$rmssdrecovery_MET_cent) + 3*sd(dataset_rmssdbase$rmssdrecovery_MET_cent)))] <- (mean(dataset_rmssdbase$rmssdrecovery_MET_cent) + 3*sd(dataset_rmssdbase$rmssdrecovery_MET_cent))
#RMSSD Recovery
dataset_rmssdbase$rmssd_recoveryno999[which(dataset_rmssdbase$rmssd_recoveryno999 < (mean(dataset_rmssdbase$rmssd_recoveryno999) - 3*sd(dataset_rmssdbase$rmssd_recoveryno999)))] <- (mean(dataset_rmssdbase$rmssd_recoveryno999) - 3*sd(dataset_rmssdbase$rmssd_recoveryno999))
dataset_rmssdbase$rmssd_recoveryno999[which(dataset_rmssdbase$rmssd_recoveryno999 > (mean(dataset_rmssdbase$rmssd_recoveryno999) + 3*sd(dataset_rmssdbase$rmssd_recoveryno999)))] <- 82

#################################################################################################################
#Basemodell RMSSD
########################################################################################################################

#Resting, Reactivity, Interaktion und Time.
model_base <- glmer(rmssd_recoveryno999 ~ resting_rmssd_log_cent + resting_met_cent + rmssdrecovery_MET_cent + rmssd_reactivity_dichotom + resting_rmssd_log_cent*rmssd_reactivity_dichotom + conttime_cent + (1|id),data = dataset_rmssdbase, family = "poisson")
tab_model(model_base, show.se = TRUE)

################################################################################################################
#Data set covariates
########################################################################################################################

dataset_rmssdco <- dataset[ , c("id", "stressmarker", "resting_rmssd", "reactivity_rmssd", "rmssd_recovery", "conttime", "resting_met", "reactivity_met", "rmssdrecovery_MET", "resting_supine", "reactivity_supine", "rmssdrecovery_supine", "rmssdrecovery_anticipatedstress", "rmssdrecovery_ambientnoise", "mvpa", "bmi", "chronicstress", "sex", "nicotine", "caffeine", "alcohol")]
dataset_rmssdco <- dataset_rmssdco[complete.cases(dataset_rmssdco),]

#RMSSD reactivity
dataset_rmssdco$rmssd_reactivity <- (dataset_rmssdco$reactivity_rmssd - dataset_rmssdco$resting_rmssd)

#RMSSD reactivity dichotom
dataset_rmssdco$rmssd_reactivity_dichotom <- dataset_rmssdco$rmssd_reactivity
dataset_rmssdco$rmssd_reactivity_dichotom [which(dataset_rmssdco$rmssd_reactivity_dichotom < 0)] <- -1
dataset_rmssdco$rmssd_reactivity_dichotom [which(dataset_rmssdco$rmssd_reactivity_dichotom > 0)] <- 1

#RMSSD Recovery
dataset_rmssdco$rmssd_recoveryno999 <- dataset_rmssdco$rmssd_recovery
dataset_rmssdco$rmssd_recoveryno999 [which(dataset_rmssdco$rmssd_recoveryno999 > 900)] <- NA
dataset_rmssdco <- dataset_rmssdco[complete.cases(dataset_rmssdco),]

#Log Transform RMSSD Parameters
#RMSSD Resting
dataset_rmssdco$resting_rmssd_log <- log(dataset_rmssdco$resting_rmssd)

#supine position in minutes
dataset_rmssdco$resting_supine_minutes <- (dataset_rmssdco$resting_supine)*5
dataset_rmssdco$reactivity_supine_minutes <- (dataset_rmssdco$reactivity_supine)*5
dataset_rmssdco$rmssdrecovery_supine_minutes <- (dataset_rmssdco$rmssdrecovery_supine)*5

#Center independent variables
#RMSSD Resting
dataset_rmssdco$resting_rmssd_log_cent <- center(dataset_rmssdco$resting_rmssd_log, type = c("CWC"), cluster = dataset_rmssdco$id, value = NULL, as.na = NULL,
                                                   check = TRUE)
#Conttime_hour
dataset_rmssdco$conttime_cent <- scale(dataset_rmssdco$conttime, scale = FALSE)
#Resting_MET
dataset_rmssdco$resting_met_cent <- scale(dataset_rmssdco$resting_met, scale = FALSE)
#Reactivity_MET
dataset_rmssdco$reactivity_met_cent <- scale(dataset_rmssdco$reactivity_met, scale = FALSE)
#Rmssdrecovery_MET
dataset_rmssdco$rmssdrecovery_MET_cent <- scale(dataset_rmssdco$rmssdrecovery_MET, scale = FALSE)
#MVPA
dataset_rmssdco$mvpa_cent <- scale(dataset_rmssdco$mvpa, scale = FALSE)
#BMI
dataset_rmssdco$bmi_cent <- scale(dataset_rmssdco$bmi, scale = FALSE)
#Chronicstress
dataset_rmssdco$chronicstress_cent <- scale(dataset_rmssdco$chronicstress, scale = TRUE)
#rmssdrecovery_supine
dataset_rmssdco$rmssdrecovery_supine_cent <- scale(dataset_rmssdco$rmssdrecovery_supine, scale = FALSE)

#Rescale variables
dataset_rmssdco$conttime_cent <- dataset_rmssdco$conttime_cent/24
dataset_rmssdco$mvpa_cent <- dataset_rmssdco$mvpa_cent/60

#outliers
#Resting Supine
dataset_rmssdco$resting_supine_minutes[which(dataset_rmssdco$resting_supine_minutes < (mean(dataset_rmssdco$resting_supine_minutes) - 3*sd(dataset_rmssdco$resting_supine_minutes)))] <- (mean(dataset_rmssdco$resting_supine_minutes) - 3*sd(dataset_rmssdco$resting_supine_minutes))
dataset_rmssdco$resting_supine_minutes[which(dataset_rmssdco$resting_supine_minutes > (mean(dataset_rmssdco$resting_supine_minutes) + 3*sd(dataset_rmssdco$resting_supine_minutes)))] <- (mean(dataset_rmssdco$resting_supine_minutes) + 3*sd(dataset_rmssdco$resting_supine_minutes))
#Reactivity Supine
dataset_rmssdco$reactivity_supine_minutes[which(dataset_rmssdco$reactivity_supine_minutes < (mean(dataset_rmssdco$reactivity_supine_minutes) - 3*sd(dataset_rmssdco$reactivity_supine_minutes)))] <- (mean(dataset_rmssdco$reactivity_supine_minutes) - 3*sd(dataset_rmssdco$reactivity_supine_minutes))
dataset_rmssdco$reactivity_supine_minutes[which(dataset_rmssdco$reactivity_supine_minutes > (mean(dataset_rmssdco$reactivity_supine_minutes) + 3*sd(dataset_rmssdco$reactivity_supine_minutes)))] <- (mean(dataset_rmssdco$reactivity_supine_minutes) + 3*sd(dataset_rmssdco$reactivity_supine_minutes))
#RMSSD Resting
dataset_rmssdco$resting_rmssd_log_cent[which(dataset_rmssdco$resting_rmssd_log_cent < (mean(dataset_rmssdco$resting_rmssd_log_cent) - 3*sd(dataset_rmssdco$resting_rmssd_log_cent)))] <- (mean(dataset_rmssdco$resting_rmssd_log_cent) - 3*sd(dataset_rmssdco$resting_rmssd_log_cent))
dataset_rmssdco$resting_rmssd_log_cent[which(dataset_rmssdco$resting_rmssd_log_cent > (mean(dataset_rmssdco$resting_rmssd_log_cent) + 3*sd(dataset_rmssdco$resting_rmssd_log_cent)))] <- (mean(dataset_rmssdco$resting_rmssd_log_cent) + 3*sd(dataset_rmssdco$resting_rmssd_log_cent))
#RMSSD Recovery
dataset_rmssdco$rmssd_recoveryno999[which(dataset_rmssdco$rmssd_recoveryno999 < (mean(dataset_rmssdco$rmssd_recoveryno999) - 3*sd(dataset_rmssdco$rmssd_recoveryno999)))] <- (mean(dataset_rmssdco$rmssd_recoveryno999) - 3*sd(dataset_rmssdco$rmssd_recoveryno999))
dataset_rmssdco$rmssd_recoveryno999[which(dataset_rmssdco$rmssd_recoveryno999 > (mean(dataset_rmssdco$rmssd_recoveryno999) + 3*sd(dataset_rmssdco$rmssd_recoveryno999)))] <- (mean(dataset_rmssdco$rmssd_recoveryno999) + 3*sd(dataset_rmssdco$rmssd_recoveryno999))
#Conttime
dataset_rmssdco$conttime_cent[which(dataset_rmssdco$conttime_cent < (mean(dataset_rmssdco$conttime_cent) - 3*sd(dataset_rmssdco$conttime_cent)))] <- (mean(dataset_rmssdco$conttime_cent) - 3*sd(dataset_rmssdco$conttime_cent))
dataset_rmssdco$conttime_cent[which(dataset_rmssdco$conttime_cent > (mean(dataset_rmssdco$conttime_cent) + 3*sd(dataset_rmssdco$conttime_cent)))] <- (mean(dataset_rmssdco$conttime_cent) + 3*sd(dataset_rmssdco$conttime_cent))
#Resting_MET
dataset_rmssdco$resting_met_cent[which(dataset_rmssdco$resting_met_cent < (mean(dataset_rmssdco$resting_met_cent) - 3*sd(dataset_rmssdco$resting_met_cent)))] <- (mean(dataset_rmssdco$resting_met_cent) - 3*sd(dataset_rmssdco$resting_met_cent))
dataset_rmssdco$resting_met_cent[which(dataset_rmssdco$resting_met_cent > (mean(dataset_rmssdco$resting_met_cent) + 3*sd(dataset_rmssdco$resting_met_cent)))] <- (mean(dataset_rmssdco$resting_met_cent) + 3*sd(dataset_rmssdco$resting_met_cent))
#Reactivity_MET
dataset_rmssdco$reactivity_met_cent[which(dataset_rmssdco$reactivity_met_cent < (mean(dataset_rmssdco$reactivity_met_cent) - 3*sd(dataset_rmssdco$reactivity_met_cent)))] <- (mean(dataset_rmssdco$reactivity_met_cent) - 3*sd(dataset_rmssdco$reactivity_met_cent))
dataset_rmssdco$reactivity_met_cent[which(dataset_rmssdco$reactivity_met_cent > (mean(dataset_rmssdco$reactivity_met_cent) + 3*sd(dataset_rmssdco$reactivity_met_cent)))] <- (mean(dataset_rmssdco$reactivity_met_cent) + 3*sd(dataset_rmssdco$reactivity_met_cent))
#Rmssdrecovery_MET
dataset_rmssdco$rmssdrecovery_MET_cent[which(dataset_rmssdco$rmssdrecovery_MET_cent < (mean(dataset_rmssdco$rmssdrecovery_MET_cent) - 3*sd(dataset_rmssdco$rmssdrecovery_MET_cent)))] <- (mean(dataset_rmssdco$rmssdrecovery_MET_cent) - 3*sd(dataset_rmssdco$rmssdrecovery_MET_cent))
dataset_rmssdco$rmssdrecovery_MET_cent[which(dataset_rmssdco$rmssdrecovery_MET_cent > (mean(dataset_rmssdco$rmssdrecovery_MET_cent) + 3*sd(dataset_rmssdco$rmssdrecovery_MET_cent)))] <- (mean(dataset_rmssdco$rmssdrecovery_MET_cent) + 3*sd(dataset_rmssdco$rmssdrecovery_MET_cent))
#MVPA
dataset_rmssdco$reactivity_supine_cent[which(dataset_rmssdco$reactivity_supine_cent < (mean(dataset_rmssdco$reactivity_supine_cent) - 3*sd(dataset_rmssdco$reactivity_supine_cent)))] <- (mean(dataset_rmssdco$reactivity_supine_cent) - 3*sd(dataset_rmssdco$reactivity_supine_cent))
dataset_rmssdco$reactivity_supine_cent[which(dataset_rmssdco$reactivity_supine_cent > (mean(dataset_rmssdco$reactivity_supine_cent) + 3*sd(dataset_rmssdco$reactivity_supine_cent)))] <- (mean(dataset_rmssdco$reactivity_supine_cent) + 3*sd(dataset_rmssdco$reactivity_supine_cent))
#BMI
dataset_rmssdco$bmi_cent[which(dataset_rmssdco$bmi_cent < (mean(dataset_rmssdco$bmi_cent) - 3*sd(dataset_rmssdco$bmi_cent)))] <- (mean(dataset_rmssdco$bmi_cent) - 3*sd(dataset_rmssdco$bmi_cent))
dataset_rmssdco$bmi_cent[which(dataset_rmssdco$bmi_cent > (mean(dataset_rmssdco$bmi_cent) + 3*sd(dataset_rmssdco$bmi_cent)))] <- (mean(dataset_rmssdco$bmi_cent) + 3*sd(dataset_rmssdco$bmi_cent))
#Chronicstress
dataset_rmssdco$chronicstress_cent[which(dataset_rmssdco$chronicstress_cent < (mean(dataset_rmssdco$chronicstress_cent) - 3*sd(dataset_rmssdco$chronicstress_cent)))] <- (mean(dataset_rmssdco$chronicstress_cent) - 3*sd(dataset_rmssdco$chronicstress_cent))
dataset_rmssdco$chronicstress_cent[which(dataset_rmssdco$chronicstress_cent > (mean(dataset_rmssdco$chronicstress_cent) + 3*sd(dataset_rmssdco$chronicstress_cent)))] <- (mean(dataset_rmssdco$chronicstress_cent) + 3*sd(dataset_rmssdco$chronicstress_cent))
#Resting Supine
dataset_rmssdco$resting_supine[which(dataset_rmssdco$resting_supine < (mean(dataset_rmssdco$resting_supine) - 3*sd(dataset_rmssdco$resting_supine)))] <- (mean(dataset_rmssdco$resting_supine) - 3*sd(dataset_rmssdco$resting_supine))
dataset_rmssdco$resting_supine[which(dataset_rmssdco$resting_supine > (mean(dataset_rmssdco$resting_supine) + 3*sd(dataset_rmssdco$resting_supine)))] <- (mean(dataset_rmssdco$resting_supine) + 3*sd(dataset_rmssdco$resting_supine))
#Reactivity Supine
dataset_rmssdco$reactivity_supine[which(dataset_rmssdco$reactivity_supine < (mean(dataset_rmssdco$reactivity_supine) - 3*sd(dataset_rmssdco$reactivity_supine)))] <- (mean(dataset_rmssdco$reactivity_supine) - 3*sd(dataset_rmssdco$reactivity_supine))
dataset_rmssdco$reactivity_supine[which(dataset_rmssdco$reactivity_supine > (mean(dataset_rmssdco$reactivity_supine) + 3*sd(dataset_rmssdco$reactivity_supine)))] <- (mean(dataset_rmssdco$reactivity_supine) + 3*sd(dataset_rmssdco$reactivity_supine))
#rmssdrecovery_supine
dataset_rmssdco$rmssdrecovery_supine[which(dataset_rmssdco$rmssdrecovery_supine < (mean(dataset_rmssdco$rmssdrecovery_supine) - 3*sd(dataset_rmssdco$rmssdrecovery_supine)))] <- (mean(dataset_rmssdco$rmssdrecovery_supine) - 3*sd(dataset_rmssdco$rmssdrecovery_supine))
dataset_rmssdco$rmssdrecovery_supine[which(dataset_rmssdco$rmssdrecovery_supine > (mean(dataset_rmssdco$rmssdrecovery_supine) + 3*sd(dataset_rmssdco$rmssdrecovery_supine)))] <- (mean(dataset_rmssdco$rmssdrecovery_supine) + 3*sd(dataset_rmssdco$rmssdrecovery_supine))

#################################################################################################################
#Model covariates
########################################################################################################################

model_covariates <- glmer(rmssd_recoveryno999 ~ resting_rmssd_log_cent + rmssd_reactivity_dichotom + resting_rmssd_log_cent*rmssd_reactivity_dichotom + conttime_cent + mvpa_cent + bmi_cent + chronicstress_cent + resting_met_cent + rmssdrecovery_MET_cent + rmssdrecovery_ambientnoise + (1|id),data = dataset_rmssdco, na.action = 'na.omit', family = "poisson")
tab_model(model_covariates, show.se = TRUE)

#################################################################################################################
#Figures
########################################################################################################################

dataset_rmssdbase$pred_mslopes <- predict(model_base)

subsetlowreact <- subset(dataset_rmssdbase, rmssd_reactivity_dichotom==0)
subsethighreact <- subset(dataset_rmssdbase, rmssd_reactivity_dichotom>0)

plot_reactivitydown <- ggplot(data=subsetlowreact, aes(x=resting_rmssd_log_cent, y=pred_mslopes, group=factor(id), colour="gray"), legend=FALSE) +
  geom_smooth(method=lm, se=FALSE, fullrange=TRUE, lty=1, size=.5, color="gray40") +
  geom_smooth(aes(group=1), method=lm, se=FALSE, fullrange=FALSE, lty=1, size=2, color="red") +
  xlab("Resting lnRMSSD") + ylab("Predicted RMSSD Recovery (min)") +
  theme_classic() +
  theme(axis.title=element_text(size=18),
        axis.text=element_text(size=18),
        plot.title=element_text(size=18, hjust=.5)) 
plot_reactivitydown <- plot_reactivitydown + coord_cartesian(xlim = c(-1, 1), ylim = c(0, 6))

plot_reactivityup <- ggplot(data=subsethighreact, aes(x=resting_rmssd_log_cent, y=pred_mslopes, group=factor(id), colour="gray"), legend=FALSE) +
  geom_smooth(method=lm, se=FALSE, fullrange=TRUE, lty=1, size=.5, color="gray40") +
  geom_smooth(aes(group=1), method=lm, se=FALSE, fullrange=FALSE, lty=1, size=2, color="blue") +
  xlab("Resting lnRMSSD") + ylab("Predicted RMSSD Recovery (min)") +
  theme_classic() +
  theme(axis.title=element_text(size=18),
        axis.text=element_text(size=18),
        plot.title=element_text(size=18, hjust=.5)) 
plot_reactivityup <- plot_reactivityup + coord_cartesian(xlim = c(-1, 1), ylim = c(0, 6))

figure <- ggarrange(plot_reactivitydown, plot_reactivityup,
                    labels = c("(A)", "(B)"),
                    ncol = 2, nrow = 1)

show(figure)
